DALS005-统计推断(Inference)04-Monte Carlo

发表于 | 分类于
| 字数统计: 2,212 | 阅读时长 ≈ 9

前言

这篇笔记是《Data Analysis for the Life Sciences》的第2章：统计推断(Inference)的第4部分。这一部分的主要内容涉及Monte Carlo模拟。

随机数

计算机可以生成伪随机数(pseudo-random numbers)，使用这些伪随机数可以用于模拟来自真实世界的随机变量。这可以让我们利用计算机来，而不是理论分析或推导来研究一些随机变量的特异。使用这种理念的一个非常有用的地方就在于，我们可以创建一些模拟数据来测试来验证我们的一些想法，而不必在实验室中做一些实验。

模拟数据也可以用于检测一些理论或分析结果的可靠性。此外，在统计学中，我们使用的很多结果是渐时性的：只当样本达到无穷大时它们才有可能成立。在实际中，我们无法收集所有的样本来验证理论与实际的关系如何。有时候我们就需要一些模拟数据。

什么是Monte Carl模拟

Monte Carl中文译为蒙特卡罗，蒙特卡罗是一种统计学方法，它使用大量的模拟数据来逼近真相，我也不太好解释，我们通过几个例子来说明一下，这一部分书本上没有，而是参考了Count Bayesie的博客，里面有不少介绍统计学的知识。

案例之一：积分

现在我们使用了一个均值为1，标准差为10的标准正态分布来讲解，下面这张图就是这个正态分布，蓝色部分是3到6构成的面积（其实就是积分），如下所示：

此时，如果我们学了微积分，就知道很容易计算出这个面积，但是，在这个案例中，我们使用Monte Carlo模拟来计算这个面积，于是，我们就能简单地使用R来写几行代码，从这个分布中采样10万次，看看有多少个值介于3到6之间，如下所示：

1
2
3
4
runs <- 100000
sims <- rnorm(runs,mean=1,sd=10)
mc.integral <- sum(sims >= 3 & sims <= 6)/runs
mc.integral

结果如下所示：

1
2
> mc.integral
[1] 0.11288

现在我们直接使用R中的累积分布函数pnorm()来计算一下，如下所示：

1
pnorm(6, mean = 1, sd = 10) -pnorm(3, mean = 1, sd = 10)

计算结果如下所示：

1
2
> pnorm(6, mean = 1, sd = 10) -pnorm(3, mean = 1, sd = 10)
[1] 0.1122028

从这两次的计算结果来看，Monte Carlo模拟出来的结果与我们实际计算的结果非常接近，可以说是完全相等了。

案例之二：估计二项分布

我们再来看一个案例，我们来抛一个硬币，抛10次，现在我们想知道，至少有3次正面朝上概论，这其实就是一个二项分布的问题，用传统的二项分布很容易计算。但是，这里我们假设我们不知道二项分布，我们就用Monte Carlo模拟来计算一下。现在我们使用0表示反面，1表示下面，然后模拟10次抛硬币的情况，我们模拟10万次，来看一下结果如何：

1
2
3
4
5
6
7
8
runs <- 100000
one.trial <- function(){
  sum(sample(c(0,1),10,replace=T)) > 3
}
mc.binom <- sum(replicate(runs,one.trial()))/runs
mc.binom

结果如下所示：

1
2
> mc.binom
[1] 0.82844

现在我们使用R中的pbinom()函数直接计算一下，如下所示：

1
pbinom(3,10,0.5,lower.tail=FALSE)

结果如下所示：

1
2
> pbinom(3,10,0.5,lower.tail=FALSE)
[1] 0.828125

前后计算出来的结果差不多。

案例之三：计算周率

我们知道，圆的面积是$A=\pi r^2$，现在我们通过Monte Carlo模拟来计算一下 $\pi$ 值。如果我们绘制了一个正方形，这个正方形的边长为$2r$，那么此时这个正方形的面积为 $A=4r^2$ ，如下所示：

此时，我们如何得到 $\pi $ 呢？我们知道，上面这个圆的面积与正方形的面积之比就是$\pi /4$ 。那么我们就可以通过一定的计算来获知这个比值的近似值，然后乘以4就得到了 $\pi$ 的近似值。接着我们就来做一下这个近似计算。

现在，我们假设上面的这个图形在一个坐标系中，圆心就是圆点 $(0,0)$ 。我们以圆心为中心，随机抽一个点，这个点的坐标是 $(x,y)$ ，圆的半径为 $r$ ，随机样本一个点后，如果 $x^2+y^2 \leq r^2$ ，那么这个点就在圆内，否则就在圆外，然后我们做10万次这样的抽取，我们计算出圆内点的数目与总的点的数目之比，然后再乘以4，就得到了 $\pi$ 的近似值，计算结果如下所示：

1
2
3
4
5
6
7
8
9
10
runs <- 100000
#runif samples from a uniform distribution
xs <- runif(runs,min=-0.5,max=0.5)
ys <- runif(runs,min=-0.5,max=0.5)
in.circle <- xs^2 + ys^2 <= 0.5^2
mc.pi <- (sum(in.circle)/runs)*4
mc.pi
plot(xs,ys,pch='.',col=ifelse(in.circle,"blue","grey")
     ,xlab='',ylab='',asp=1,
     main=paste("MC Approximation of Pi =",mc.pi))

计算结果如下所示：

1
2
> mc.pi
[1] 3.14124

从计算结果来看，我们得到的 $\pi$ 近似为3.14124，还是比较接近的。

Monte Carlo模拟

现在我们再回到原书，我们使用Monte Carlo模拟来比较一下在不同样本数目中，CLT与t分布的接近程度，如下所示：

1
2
3
4
5
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist

现在构建一个函数，这个函数用于生成零假设下，样本数目n的t统计量，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/mice_pheno.csv")
dat <- read.csv(filename)
controlPopulation <- filter(dat, Sex == "F" & Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist
ttestgenerator <- function(n){
  # note that here we have a false "high fat" group where we actually
  # sample from the nonsmokers. this is because we are modeling the *null*
  cases <- sample(controlPopulation, n)
  controls <- sample(controlPopulation,n)
  tstat <- (mean(cases) - mean(controls))/
    sqrt(var(cases)/n + var(controls)/n)
  return(tstat)
}
ttests <- replicate(1000,ttestgenerator(10))
hist(ttests)

这个t分布是否非常接近于正态分布？我们可以使用qq图(quantile-quantile)来看一下，如下所示：

1
2
qqnorm(ttests)
abline(0,1)

从上图我们可以看出来，这个数据集非常接近于正态分布。如果现在我们将样本数改为3，qq图如下所示：

1
2
3
ttests <- replicate(1000, ttestgenerator(3))
qqnorm(ttests)
abline(0,1)

从上图可以看出来，数据的高分位数(large quantiles)区，也就是统计学家称的尾部(tails)，比较大，这个尾部区就是上图中直线左侧的下方数据，以及右侧直线上方的数据。在前面的部分我们提到，即使总体服从正态分布，那么小样本则是更加近似地服从t分布，下面我们来的模拟一下：

1
2
3
ps <- (seq(0,999)+0.5)/1000
qqplot(qt(ps,df=2*3-2),ttests,xlim=c(-6,6),ylim=c(-6,6))
abline(0,1)

从上面两张图可以看出来，数据更接近于t分布，而不是正态分布，不过虽然接近于t分布，但是这种接近也并非完全，这是由于原始数据也并不是完美的服从正态分布，如下所示：

1
2
qqnorm(controlPopulation)
qqline(controlPopulation)

观测值的参数化模拟(Parametric Simulations for the Observations)

上述我们模拟随机变量的方法叫Monte Carlo模拟。我们利用既有的总体数据，随机抽取了一些样本。在实际分析中，我们无法获取整个总体。但是，当我们想在实际分析过程中使用Monte Carlo模拟时，一种比较典型的做法是，假设一个参数分布(parametric distribution)，并根据此分布生成一个总体，这种手段叫参数模拟(parametric simulation)。这就意味着，我们从真实的数据（这里指的是样本的均值与标准差）中提取参数，将这些参数添加一个模型（这里指的是正态分布）。这是Monte Carlo模拟最普遍的手段。

回到小鼠的体重案例上来，我们可以利用我们现在有的知识，例如小鼠的体重通常是24g，SD为3.5g，小鼠的体重用了近服从正态分布，我们可以产生一个总体：

1
controls<- rnorm(5000, mean=24, sd=3.5)

当我们生成了这批数据后，我们可以重复再生成几次。此时，我们并不使用sample()函数来进行抽样，而是使用随机数来生成一些服从正态分布的数据，前面提到的ttestgenerator()函数更改如下：

1
2
3
4
5
6
7
ttestgenerator <- function(n, mean=24, sd=3.5){
  cases <- rnorm(n, mean,sd)
  controls <- rnorm(n, mean, sd)
  tstat <- (mean(cases)-mean(controls))/
    sqrt(var(cases)/n + var(controls)/n)
  return(tstat)
}

练习

P88